function result = mie_abcd(m, x)

  % Computes a matrix of Mie coefficients, a_n, b_n, c_n, d_n,
  % of orders n=1 to nmax, complex refractive index m=m'+im",
  % and size parameter x=k0*a, where k0= wave number
  % in the ambient medium, a=sphere radius;
  % p. 100, 477 in Bohren and Huffman (1983) BEWI:TDD122
  % C. Mätzler, June 2002

  nmax = round(2 + x + 4 * x.^(1/3))
  n = (1:nmax)
  nu = (n + 0.5)
  z = m .* x
  m2 = m .* m
  sqx = sqrt(0.5 * pi ./ x)
  sqz = sqrt(0.5 * pi ./ z)
  bx = besselj(nu, x) .* sqx
  bz = besselj(nu, z) .* sqz
  yx = bessely(nu, x) .* sqx
  hx = bx + i * yx
  b1x = [sin(x) / x, bx(1:nmax - 1)]
  b1z = [sin(z) / z, bz(1:nmax - 1)]
  y1x = [-cos(x) / x, yx(1:nmax - 1)]
  h1x = b1x + i * y1x
  ax = x .* b1x - n .* bx
  az = z .* b1z - n .* bz
  ahx = x .* h1x - n .* hx

  an = (m2 .* bz .* ax - bx .* az) ./ (m2 .* bz .* ahx - hx .* az)
  bn = (bz .* ax - bx .* az) ./ (bz .* ahx - hx .* az)
  cn = (bx .* ahx - hx .* ax) ./ (bz .* ahx - hx .* az)
  dn = m .* (bx .* ahx - hx .* ax) ./ (m2 .* bz .* ahx - hx .* az)
  result = [an; bn; cn; dn]
